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Abstract 

Time integration of ODEs or time-dependent PDEs with required resolution of the 
fastest time scales of the system, can be very costly if the system exhibits multiple 
time scales of different magnitudes. If the different time scales are localised to dif- 
ferent components, corresponding to localisation in space for a PDE, efficient time 
integration thus requires that we use different time steps for different components. 

We present an overview of the multi-adaptive Galerkin methods mcG(g) and 
mdG(g) recently introduced in a series of papers by the author. In these methods, 
the time step sequence is selected individually and adaptively for each component, 
based on an a posteriori error estimate of the global error. 

The multi-adaptive methods require the solution of large systems of nonlinear 
algebraic equations which are solved using explicit-type iterative solvers (fixed point 
iteration). If the system is stiff, these iterations may fail to converge, corresponding 
to the well-known fact that standard explicit methods are inefficient for stiff systems. 
To resolve this problem, we present an adaptive strategy for explicit time integration 
of stiff ODEs, in which the explicit method is adaptively stabilised by a small number 
of small, stabilising time steps. 
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1 Introduction 



In earlier work [29,30], we have introduced the mult i- adaptive Galerkin meth- 
ods mcG(g) and mdG(g) for ODEs of the type 

' u(t) = f(u(t),t), te(0,T], 
^ u(0) = u , 

where u : [0, T] — > R N is the solution to be computed, u G R N a given 
initial condition, T > a given final time, and / : R N x (0,T] ->■ R N a 
given function that is Lipschitz-continuous in u and bounded. We use the 
term multi-adaptivity to describe methods with individual time-stepping for 
the different components Ui(t) of the solution vector u(t) = (ui(t)), including 
(i) time step length, (ii) order, and (iii) quadrature, all chosen adaptively in 
a computational feed-back process. 

Surprisingly, individual time-stepping for ODEs has received little attention 
in the large literature on numerical methods for ODEs, see e.g. [4,24,25,2,34]. 
For specific applications, such as the n-body problem, methods with individual 
time-stepping have been used, see e.g. [31,1,5], but a general methodology has 
been lacking. For time-dependent PDEs, in particular for conservation laws 
of the type u + f(u) x = 0, attempts have been made to construct methods 
with individual (locally varying in space) time steps. Flaherty et al. [22] have 
constructed a method based on the discontinuous Galerkin method combined 
with local explicit Euler time-stepping. A similar approach is taken in [6] where 
a method based on the original work by Osher and Sanders [32] is presented 
for conservation laws in one and two space dimensions. Typically the time 
steps used are based on local CFL conditions rather than error estimates for 
the global error and the methods are low order in time (meaning < 2). We 
believe that our work on multi-adaptive Galerkin methods (including error 
estimation and arbitrary order methods) presents a general methodology to 
individual time-stepping, which will result in efficient integrators both for 
ODEs and time-dependent PDEs. 

The mult i- adaptive methods are developed within the general framework of 
adaptive Galerkin methods based on piecewise polynomial approximation (fi- 
nite element methods) for differential equations, including the continuous and 
discontinuous Galerkin methods cG(q) and dG(q) which we extend to their 
multi-adaptive analogues mcG(g) and mdG(g). Earlier work on adaptive error 
control for the cG(q) and dG(q) methods include [7,17,27,19,18,21]. See also 
[10,11,9,12,13,14], and [8] or [20] in particular for an overview of adaptive error 
control based on duality techniques. The approach to error analysis and adap- 
tivity presented in these references naturally carries over to the mult i- adaptive 
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methods. 



1.1 The stiffness problem 

The classical wisdom developed in the 1950s regarding stiff ODEs is that ef- 
ficient integration requires implicit (A-stable) methods, at least outside tran- 
sients where the time steps may be chosen large from accuracy point of view. 
Using an explicit method (with a bounded stability region) the time steps 
have to be small at all times for stability reasons, in particular outside tran- 
sients, and the advantage of a low cost per time step for the explicit method 
is counter-balanced by the necessity of taking a large number of small time 
steps. As a result, the overall efficiency of an explicit method for a stiff ODE 
is small. 

We encounter the same problem when we try to use explicit fixed point it- 
eration to solve the discrete equations given by the multi-adaptive Galerkin 
methods mcG(g) and mdG(g). However, it turns out that if a sequence of 
large (unstable) time steps are accompanied by a suitable (small) number of 
small time steps, a stiff system can be stabilised to allow integration with 
an effective time step much larger than the largest stable time step given by 
classical stability analysis. This idea of stabilising a stiff system using the in- 
herent damping property of the stiff system itself was first developed in an 
automatic and adaptive setting in [16], and will be further explored in the 
full mult i- adaptive setting. A similar approach is taken in recent independent 
work by Gear and Kevrekidis [3] . The relation to Runge-Kutta methods based 
on Chebyshev polynomials discussed by Verwer in [26] should also be noted. 



1.2 Notation 

The following notation is used in the discussion of the multi-adaptive Galerkin 
methods below: Each component Ui(t), i = 1,...,N, of the approximate 
m(c/d)G(g) solution U(t) of (1) is a piecewise polynomial on a partition 
of (0, T] into Mi subintervals. Subinterval j for component % is denoted by 
hj = (ti,j-i,Uj], and the length of the subinterval is given by the local time 
step kij = tij — Uj-i. This is illustrated in Figure 1. On each subinterval fj, 
Ui\j tj is a polynomial of degree q^ and we refer to (1^, U^i^) as an element. 

Furthermore, we shall assume that the interval (0, T] is partitioned into blocks 
between certain synchronised time levels = To < T\ < . . . < Tm = T. We 
refer to the set of intervals T n between two synchronised time levels T n _i and 
T n as a time slab: T n = {fj : T n _x < Uj-i < t^ < T n }, and we denote the 
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length of a time slab by K n = T n — T n _\. The partition consisting of the entire 
collection of intervals is denoted by T = U7^. 
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Fig. 1. Individual partitions of the interval (0, T] for different components. Ele- 
ments between common synchronised time levels are organised in time slabs. In this 
example, we have N = 6 and M = 4. 



1.3 Outline 



The outline of the paper is as follows: In Section 2, we formulate the multi- 
adaptive Galerkin methods mcG(g) and mdG(g). In Section 3, we discuss error 
control and adaptivity. In particular, we show how to choose the individual 
time steps based on an a posteriori error estimate for the global error. In 
Section 4, we give a quick overview of an iterative method (based on fixed 
point iteration) for the system of nonlinear discrete equations that needs to 
be solved on each time slab, and in Section 5 we describe a technique that 
can be used to stabilise the explicit fixed point iterations for stiff problems. 
Finally, in Section 6, we present a number of numerical examples chosen to 
illustrate both the potential of multi-adaptivity and the use of explicit fixed 
point iteration (or explicit time-stepping) for stiff problems. 



4 



2 Multi-Adaptive Galerkin 



2.1 Multi- adaptive continuous Galerkin, mcG(g) 

To give the definition of the mcG(g) method, we define the trial space V and 
the test space W as 

V = {ve [C([0, T])] N : Vi \ hj G V^(Iij), j = 1, . . • , Mi, i — 1,. . ., N}, 
W = {v:v i \ Iij eV**- 1 (I ij ), j = l,...,Mi, i = l,...,N}, 

where V q (I) denotes the linear space of polynomials of degree q > on the 
interval /. In other words, V is the space of continuous piecewise polynomials 
of degree q = qiit) = q^, t G 1^ on the partition T, and W is the space of 
(in general discontinuous) piecewise polynomials of degree q — 1 on the same 
partition. 

We define the mcG(g) method for (1) as follows: Find U G V with U(0) = uo, 
such that 

T T 

j(U,v)alt = J(f(U,-),v)dt VveW, (3) 





where (•, •) denotes the standard inner product in WL N . If now for each local 
interval we take v n = when n ^ i and Vi(t) = when t £ Iij, we can 
rewrite the global problem (3) as a number of successive local problems for 
each component: For i = 1, . . . , N, j = 1, . . . , Mj, find U^i^ G V qij (Iij) with 
Ui(tij-i) given from the previous time interval, such that 

J UiV dt = J fi(U, > dt \/v G V" '(/;,). (4) 



We define the residual R of the approximate solution U to be R(U, t) — U (t) — 
f(U(t),t). In terms of the residual, we can rewrite (4) as J t Ri(U, -)v dt = 
for all v G V qii ~ l {Iij), i.e., the residual is orthogonal to the test space on every 
local interval. We refer to this as the Galerkin orthogonality of the mcG(g) 
method. 
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2.2 Multi- adaptive discontinuous Galerkin, mdG(g) 



For the mdG(g) method, we define the trial and test spaces by 

V = W = {v : v t \ hj G V"Hl tJ ).j = 1, • • • , Mi, i — 1, . . . , N}, (5) 



i.e., both trial and test functions are (in general discontinuous) piecewise poly- 
nomials of degree q = q^t) = qij, t G Iy on the partition T. 

We define the mdG(g) method for (1) as follows, similar to the definition of 
the continuous method: Find U G V with U(0~) = u , such that 



N Mi 

EE 

i=i j=i 



[u^j.m^-i) + J dt 
i, 



(f(U,-),v)dt VveW. (6) 



where [•] denotes the jump, i.e., [v]ij = v(tfj) — vfj:^). 

The mdG(g) method in local form, corresponding to (4), reads: For i = 
1, . . . , N, j = 1, . . . , M h find Ui\ hj G V qi i (Iij), such that 



[^L-i^j-i) + Ju t vdt = J f t (U, > dt Mv G V^ihj), 



(7) 



where the initial condition is specified for i = 1,...,N, by Ui(0~) = Mj(0). 
In the same way as for the continuous method, we define the residual R of 
the approximate solution U to be R{U,t) = U(t) — f(U(t),t), defined on the 
inner of every local interval i^, and rewrite (7) in the form [Ui\ij-iv(ti t j-i) + 
jj..Ri(U,-)v dt = for all v G V q ^(Iij). We refer to this as the Galerkin 
orthogonality of the mdG(g) method. 



3 Error Control and Adaptivity 

Our goal is to compute an approximation U(T) of the exact solution u{T) 
of (1) at final time T within a given tolerance TOL > 0, using a minimal 
amount of computational work. This goal includes an aspect of reliability (the 
error should be less than the tolerance) and an aspect of efficiency (minimal 
computational work). To measure the error we choose a norm, such as the 
Euclidean norm || • || on M. N , or more generally some other quantity of interest. 
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We discuss below both a priori and a posteriori error estimates for the multi- 
adaptive Galerkin methods, and the application of the a posteriori error esti- 
mates in multi-adaptive time-stepping. 



3.1 A priori error estimates 



Standard (duality-based) a priori error estimates show that the order for the 
ordinary Galerkin methods cG(g) and dG(g) is 2q and 2q + 1, respectively. 
A generalisation of these estimates to the multi-adaptive methods gives the 
same result. The mult i- adaptive continuous Galerkin method mcG(g) is thus 
of order 2q, and the mult i- adaptive discontinuous Galerkin method mdG(g) is 
of order 2q + 1. 



3.2 A posteriori error estimates 



A posteriori error analysis in the general framework of [8] relies on the concept 
of the dual problem. The dual problem of the initial value problem (1) is the 
linearised backward problem given by 



-<}> = J*(u,U,-)(f> on[0,T), 
d>(T)=e(T)/\\e(T)l 



where the Jacobian J is given by J(u, U, •) = Jq 1 |^ (su + (1 — s)U, ■) ds and 
* denotes the transpose. We use the dual problem to represent the error in 
terms of the dual solution and the residual R. For the mcG(g) method the 
representation formula is given by 



T 



;(T)|| = J(RA)dt, (9) 



and for the mdG(g) method, we obtain 

N Mi 

l|e(T)|| = EE Pih-iUkj-i) + / MU,-)cf>i dt. (10) 
i=ij=i I 

1 % i 
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Using the Galerkin orthogonalities together with special interpolation esti- 
mates (see [29]), we obtain a posteriori error estimates of the form 

N 

\\e(T)\\<Y,Sl qA m^{CkTn}, (11) 
i=i ^ U ' J ' 



for the mcG(g) method, and 

N 

\\e(T)\\<^St +1] m^{Ckr +l ri }, (12) 



i=l 



for the mdG(g) method, where C is an interpolation constant, r,i is a local 
measure of the residual, and the individual stability factors Si are given by 
S\ 9 ^ = Jq \<f)\ q ^\ dt. Typically, the stability factors are of moderate size for a 
stiff problem (and of unit size for a parabolic problem), which means that ac- 
curate computation is possible over long time intervals. Note that the Lipschitz 
constant, which is large for a stiff problem, is not present in these estimates. 

The analysis can be extended to include also computational errors, arising 
from solving the discrete equations using an iterative method, and quadrature 
errors, arising from evaluating the integrals in (4) and (7) using quadrature. 



3. 3 Adaptivity 



To achieve the goals stated at the beginning of this section, the adaptive 
algorithm chooses individual time steps for the different components based on 
the a posteriori error estimates. Using for example a standard PID regulator 
from control theory, we choose the individual time steps for each component 
to satisfy 

S.Ck-j'v,, = TOL/iV, (13) 



or, taking the logarithm with Q = log(TOL/(7VS;C)), 

Pij log kij + log nj = Q, (14) 

with maximal time steps {%}, following work by Soderlind and coworkers 
[23,35]. Here, p^ = qij for the mcG(g) method and p^ = qij + 1 for the 
mdG(g) method. 

To solve the dual problem (8), which is needed to compute the stability factors, 
it would seem that we need to know the error e(T), since this is used as an 
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initial value for the dual problem. However, we know from experience that 
the stability factors are quite insensitive to the choice of initial data for the 
dual problem. (A motivation of this for parabolic problems is given in [15].) 
Thus in practice, a random (and normalised) value is chosen as initial data for 
the dual problem. Another approach is to take the initial value for the dual 
problem to be (f>(T) = (0, . . . , 0, 1, 0, . . . , 0), i.e., a vector of zeros except for 
a single component which is of size one. This gives an estimate for a chosen 
component of the error. By other choices of data for the dual problem, other 
functional of the error can be controlled. In either case, the stability factors 
are computed using quadrature from the computed dual solution. 

The adaptive algorithm can be expressed as follows: Given a tolerance TOL > 
0, make a preliminary estimate for the stability factors and then 

(i) Solve the primal problem with time steps based on (13); 

(ii) Solve the dual problem and compute the stability factors; 

(iii) Compute an error bound E based on (9) or (10); 

(iv) If E < TOL then stop; if not go back to (i). 

Note that we use the error representations (9) and (10) to obtain sharp error 
estimates. On the other hand, the error estimates (11) and (12) are used to 
determine the adaptive time step sequences. 

To limit the computational work, it is desirable that only a few iterations in 
the adaptive algorithm are needed. In the simplest case, the error estimate will 
stay below the given tolerance on the first attempt. Otherwise, the algorithm 
will try to get below the tolerance the second time. It is also possible to limit 
the number of times the dual problem is solved. It should also be noted that 
to obtain an error estimate at a time t = i , different from the final time T, 
the dual problem has to be solved backwards also from time i. This may be 
necessary in some cases if the stability factors do not grow monotonically as 
functions of the final time T, but for many problems the stability factors grow 
with T, indicating accumulation of errors. 

Our experience is that automatic computation based on this adaptive strategy 
is both reliable (the error estimates are quite close to the actual error) and 
efficient (the additional cost for solving the dual problem is quite small). See 
[28] for a discussion on this topic. 



4 Iterative Methods for the Nonlinear System 

The nonlinear discrete algebraic equations given by the mcG(g) and mdG(g) 
methods presented in Section 2 (including numerical quadrature) to be solved 
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on every local interval 1^ take the form 

lij 

n=0 

for m — 0, . . . , g^j, where {£ijm}m=o are the degrees of freedom to be deter- 
mined for component Ui(t) on the interval 1^, {«Jmn}* 3 =0n=0 are weights, re- 
maps lij to (0, 1]: Tij(t) = (t—tij-i)/(tij—ti t j-i), and {s!?' 3 '}*i are quadrature 
points defined on [0, 1]. 

The strategy we use to solve the discrete equations (15) is by direct fixed 
point iteration, possibly in combination with a simplified Newton's method. 
To evolve the system, we need to collect the degrees of freedom for different 
components between two time levels and solve the discrete equations for these 
degrees of freedom. We refer to such a collection of elements between two time 
levels as a time slab (see Figure 1). New time slabs are formed as we evolve 
the system starting at time t = 0, in the same way as new time intervals are 
formed in a standard solver which uses the same time steps for all components. 
On each time slab, we thus compute the degrees of freedom {£ij m }m=o f° r eacn 
element within the time slab using (15), and repeat the iterations until the 
computational error is below a given tolerance for the computational error. 
The iterations are carried out in order, starting at the element closest to time 
t = and continuing until we reach the last element within the time slab. 
This is illustrated in Figure 2. 




Tn—l T n 

Fig. 2. Multi-adaptive time-stepping within a time slab for a system with two com- 
ponents. 

The motivation for using direct fixed point iteration, rather than using a full 
Newton's method, is that we want to avoid forming the Jacobian (which may 
be very large, since the nonlinear system to be solved is for the entire time 
slab) and also avoid solving the linearised system. Instead, using the strategy 
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for adaptive damping of the fixed point iterations as described in the next 
section, the linear algebra is built into the adaptive solver. Since often only a 
small number of fixed point iterations are needed, typically only two or three 
iterations, we believe this to be an efficient approach. 



5 Stiff Problems 

As discussed in the previous section, the nonlinear discrete equations given by 
the (implicit) mult i- adaptive Galerkin methods are solved using fixed point 
iteration on each time slab. For stiff problems these iterations may fail to 
converge. We now discuss a simple way to stabilise a stiff system, in order to 
make the explicit fixed point iterations convergent. 

For simplicity, we assume that the time step sequence, ki, &2, • • • , &Ar, is the 
same for all components. 



5.1 The test equation 

To demonstrate the main idea, we consider the stabilisation of the explicit 
Euler method applied to the simple test equation: 

u(t) + Xu(t) =0 for t > 0, 

(16) 

U{0) = M , 

where A > and uq is a given initial condition. The solution is given by 
u(t) = exp(— Xt)uo- 

The explicit Euler method for the test equation reads 

jjn = V n-1 _ knXU n-l = ^ _ fc^J/n-l. 

This method is conditionally stable, with stability guaranteed if k n \ < 2. If A 
is large, this is too restrictive outside transients. 

Now, let K be a large time step satisfying K\ > 2 and let k a small time step 
chosen so that k\ < 2. Consider the method 

U n = (1 - A;A) m (l -K\)U n -\ (17) 
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corresponding to one explicit Euler step with large time step K and m ex- 
plicit Euler steps with small time steps k, where m is a positive integer to be 
determined. Altogether this corresponds to a time step of size k n = K + mk. 
For the overall method to be stable, we require that |1 — kX\ m (KX — 1) < 1, 
that is 

m > y«-i) iog(jrA) - 

— log 1 1 — k\\ c 



if KX ^> 1 and c = k\ is of moderate size, say c = 1/2. 

We conclude that m will be quite small and hence the small time steps will 
be used only in a small fraction of the total time interval, giving a large 
effective time step. To see this, define the cost as a = £ (l/K,l/k), i.e., 

the number of time steps per unit interval. Classical stability analysis gives 
a — 1/k — A/2 with a maximum time step k — 2/A. Using (18) we instead 
find 



for KX ^> 1. The cost is thus decreased by the cost reduction factor 

2\og(KX) \og(KX) 
cKX KX ' 



which can be quite significant for large values of KX. 



5.2 The general non-linear problem 



For the general nonlinear problem (1), the gain is determined by the distribu- 
tion of the eigenvalues of the Jacobian, see [16]. The method of stabilising the 
system using a couple of small stabilising time steps is best suited for systems 
with a clear separation of the eigenvalues into small and large eigenvalues, but 
even for the semi-discretised heat equation (for which we have a whole range 
of eigenvalues) the gain can be substantial, as we shall see below. 



5.3 An adaptive algorithm 



In [16] we present an adaptive algorithm in which both the size of the small sta- 
bilising time steps and the number of such small time steps are automatically 
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determined. Using adaptive stabilisation, the damping is targeted precisely 
at the current unstable eigenmode, which as a consequence allows efficient 
integration also of problems with no clear separation of its eigenvalues. 



6 Numerical Examples 

The numerical examples presented in this section are divided into two cat- 
egories: examples illustrating the concept of mult i- adapt ivity and examples 
illustrating explicit time-stepping (or explicit fixed point iteration) for stiff 
problems. 



6.1 Multi-adaptivity 

The two examples presented below are taken from [30], in which further ex- 
amples are presented and discussed in more detail. 

6.1.1 A mechanical multi-scale system 

To demonstrate the potential of the multi-adaptive methods, we consider a 
dynamical system in which a small part of the system oscillates rapidly. The 
problem is to compute accurately the positions (and velocities) of the N point- 
masses attached together with springs of equal stiffness as in Figure 3. 




Fig. 3. A mechanical system consisting of TV = 5 masses attached together with 
springs. 

We choose a small time step for the smallest mass and large time steps for the 
larger masses, and measure the work for the mcG(l) method as we increase 
the number of larger masses. The work is then compared to the work required 
for the standard cG(l) method using the same (small) time step for all masses. 
As is evident in Figure 4, the work (in terms of function evaluations) increases 
linearly for the standard method, whereas for the mult i- adaptive method it 
remains practically constant. 
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Fig. 4. Error, cpu time, total number of steps, and number of function evaluations 
as function of the number of masses, for the multi-adaptive cG(l) method (dashed) 
and the standard cG(l) method (solid). 



6.1.2 Reaction- diffusion 

Next consider the following system of PDEs: 



U2 — €U 2 = Wl«2> 



(20) 



on (0,1) x (0, T] with e = 0.001, T = 100 and homogeneous Neumann 
boundary conditions at x = and x — 1, which models isothermal auto- 
catalytic reactions (see [33]): A\ + 2A 2 — > A 2 + 2A 2 . As initial conditions, 
we take U\{x, 0) = for < a; < x , Ui(a;,0) = 1 for Xq < x < 1, and 
■u 2 (x, 0) = 1 — Ui(x, 0) with x = 0.2. An initial reaction where substance A\ is 
consumed and substance A 2 is formed will then take place at x = xq, resulting 
in a decrease in the concentration u\ and an increase in the concentration u 2 . 
The reaction then propagates to the right until all of substance A\ is consumed 
and we have u\ — and u 2 = 1 in the entire domain. 

Computing the solution using the mcG(2) method, we find that the time steps 
are automatically chosen to be small only in the vicinity of the reaction front, 
see Figure 5, and during the computation the region of small time steps will 
propagate to the right at the same speed as the reaction front. 
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Fig. 5. The concentrations of the two species, U% and U2, at time t = 50 as function 
of space (above), and the corresponding time steps (below). 



6.2 Explicit time- stepping for stiff problems 



To illustrate the technique of stabilisation for stiff problems, we present below 
some examples taken from [16]. In these examples, the cost a is compared to 
the cost «o of a standard implementation of the cG(l) method in which we 
are forced to take a small time step all the time. (These small time steps are 
marked by dashed lines in the figures.) Comparison has not been made with 
an implicit method, since it would be difficult to make such a comparison fair; 
one could always argue about the choice of linear solver and preconditioner. 
However, judging by the modest restriction of the average time step size and 
the low cost of the explicit method, we believe our approach to be competitive 
also with implicit methods, although this remains to be seen. 



6.2.1 The test equation 

The first problem we try is the test equation: 

u(t) + Xu(t) =0 for t > 0, 
< W W (21) 

u(0) = u , 



on [0, 10], where we choose Uq = 1 and A = 1000. As is shown in Figure 6, the 
time step is repeatedly decreased to stabilise the stiff system, but overall the 
effective time step is large and the cost reduction factor is a/ao ~ 1/310. 
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Fig. 6. Solution and time step sequence for eq. (21), a/ao « 1/310. 



6.2.2 The test system 
For the test system, 



u(t) + Au(t) =0 for t > 0, 
W W (22) 

w(0) = u , 



on [0,10], we take A = diag(100, 1000) and -u = (1>1)- As seen in Figure 
7, most of the stabilising steps are chosen to damp out the eigenmode cor- 
responding to the largest eigenvalue, A2 = 1000, but some of the damping 
steps are targeted at the second eigenvalue, Ai = 100. The selective damping 
is handled automatically by the adaptive algorithm and the cost reduction 
factor is again significant: a/a ~ 1/104. 



6.2.3 The HIRES problem 

The so-called HIRES problem ( "High Irradiance RESponse" ) originates from 
plant physiology and is taken from the test set of ODE problems compiled by 
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Fig. 7. Solution and time step sequence for eq. (22), a/ao « 1/104. 
Lioen and de Swart [36] . The problem consists of the following eight equations: 



ill = 


-1.71«i + 0.43u 2 + 8.32u 3 + 0.0007, 


u 2 = 


1.71«i -8.75u 2 , 


U3 = 


-10.03u 3 + 0.43u 4 + 0.035m 5 , 


ii4 = 


8.32m 2 + 1.71« 3 - 1.12« 4 , 


ii 5 = 


-1.745u 5 + 0.43u 6 + 0.43m 7 , 


u 6 = 


-280.0w 6 w 8 + 0.69w 4 + 1.71u B - 0.43w 6 + 0.69m 


u 7 = 


280.0m 6 m 8 - 1.81w 7 , 


u 8 = 


-280.0m 6 m 8 + 1.81w 7 , 



(23) 



on [0,321.8122] (as specified in [36]). The initial condition is given by Uq = 
(1.0,0,0,0,0,0,0,0.0057). The cost reduction factor is now a/a ~ 1/33, see 
Figure 8. 

6.3 The heat equation 

Finally, we consider the heat equation in one dimension: 

u(x, t) - u"(x, t) = f(x, t), x e (0, 1), t > 0, 

u (0) = = 0, (24) 
u(-,t) = 0, 
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Fig. 8. Solution and time step sequence for eq. (23), a/ao « 1/33. 

where we choose f(x, t) = f(x) as an approximation of the Dirac delta function 
at x = 0.5. Discretising in space, we obtain the ODE 



u(t) + Au(t) = f, t > 0, 
u(0) = 0, 



(25) 



where A is the stiffness matrix. With a spatial resolution of h = 0.01, the 
eigenvalues of A are distributed in the interval [0,4 • 10 4 ] (see Figure 9). The 
selective damping produced by the adaptive algorithm performs well and the 
cost reduction factor is a/ao ~ 1/17. 
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